home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / bit / src / mpeg / jrevdct.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  44KB  |  1,461 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on an algorithm described in
  11.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  12.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  13.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  14.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  15.  * We use their alternate method with 12 multiplies and 32 adds.
  16.  * The advantage of this method is that no data path contains more than one
  17.  * multiplication; this allows a very simple and accurate implementation in
  18.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  19.  * 
  20.  * I've made lots of modifications to attempt to take advantage of the
  21.  * sparse nature of the DCT matrices we're getting.  Although the logic
  22.  * is cumbersome, it's straightforward and the resulting code is much
  23.  * faster.
  24.  *
  25.  * A better way to do this would be to pass in the DCT block as a sparse
  26.  * matrix, perhaps with the difference cases encoded.
  27.  */
  28.  
  29. #include <string.h>
  30. #include "video.h"
  31. #include "proto.h"
  32.  
  33. #define GLOBAL            /* a function referenced thru EXTERNs */
  34.   
  35. /* We assume that right shift corresponds to signed division by 2 with
  36.  * rounding towards minus infinity.  This is correct for typical "arithmetic
  37.  * shift" instructions that shift in copies of the sign bit.  But some
  38.  * C compilers implement >> with an unsigned shift.  For these machines you
  39.  * must define RIGHT_SHIFT_IS_UNSIGNED.
  40.  * RIGHT_SHIFT provides a proper signed right shift of an INT32 quantity.
  41.  * It is only applied with constant shift counts.  SHIFT_TEMPS must be
  42.  * included in the variables of any routine using RIGHT_SHIFT.
  43.  */
  44.   
  45. #ifdef RIGHT_SHIFT_IS_UNSIGNED
  46. #define SHIFT_TEMPS    INT32 shift_temp;
  47. #define RIGHT_SHIFT(x,shft)  \
  48.     ((shift_temp = (x)) < 0 ? \
  49.      (shift_temp >> (shft)) | ((~((INT32) 0)) << (32-(shft))) : \
  50.      (shift_temp >> (shft)))
  51. #else
  52. #define SHIFT_TEMPS
  53. #define RIGHT_SHIFT(x,shft)    ((x) >> (shft))
  54. #endif
  55.  
  56. /*
  57.  * This routine is specialized to the case DCTSIZE = 8.
  58.  */
  59.  
  60. #if DCTSIZE != 8
  61.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  62. #endif
  63.  
  64.  
  65. /*
  66.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  67.  * on each column.  Direct algorithms are also available, but they are
  68.  * much more complex and seem not to be any faster when reduced to code.
  69.  *
  70.  * The poop on this scaling stuff is as follows:
  71.  *
  72.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  73.  * larger than the true IDCT outputs.  The final outputs are therefore
  74.  * a factor of N larger than desired; since N=8 this can be cured by
  75.  * a simple right shift at the end of the algorithm.  The advantage of
  76.  * this arrangement is that we save two multiplications per 1-D IDCT,
  77.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  78.  *
  79.  * We have to do addition and subtraction of the integer inputs, which
  80.  * is no problem, and multiplication by fractional constants, which is
  81.  * a problem to do in integer arithmetic.  We multiply all the constants
  82.  * by CONST_SCALE and convert them to integer constants (thus retaining
  83.  * CONST_BITS bits of precision in the constants).  After doing a
  84.  * multiplication we have to divide the product by CONST_SCALE, with proper
  85.  * rounding, to produce the correct output.  This division can be done
  86.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  87.  * as long as possible so that partial sums can be added together with
  88.  * full fractional precision.
  89.  *
  90.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  91.  * they are represented to better-than-integral precision.  These outputs
  92.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  93.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  94.  * intermediate INT32 array would be needed.)
  95.  *
  96.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  97.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  98.  * shows that the values given below are the most effective.
  99.  */
  100.  
  101. #ifdef EIGHT_BIT_SAMPLES
  102. #define PASS1_BITS  2
  103. #else
  104. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  105. #endif
  106.  
  107. #define ONE    ((INT32) 1)
  108.  
  109. #define CONST_SCALE (ONE << CONST_BITS)
  110.  
  111. /* Convert a positive real constant to an integer scaled by CONST_SCALE.
  112.  * IMPORTANT: if your compiler doesn't do this arithmetic at compile time,
  113.  * you will pay a significant penalty in run time.  In that case, figure
  114.  * the correct integer constant values and insert them by hand.
  115.  */
  116.  
  117. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  118.  
  119. /* Descale and correctly round an INT32 value that's scaled by N bits.
  120.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  121.  * the fudge factor is correct for either sign of X.
  122.  */
  123.  
  124. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  125.  
  126. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  127.  * For 8-bit samples with the recommended scaling, all the variable
  128.  * and constant values involved are no more than 16 bits wide, so a
  129.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  130.  * this provides a useful speedup on many machines.
  131.  * There is no way to specify a 16x16->32 multiply in portable C, but
  132.  * some C compilers will do the right thing if you provide the correct
  133.  * combination of casts.
  134.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  135.  */
  136.  
  137. #ifdef EIGHT_BIT_SAMPLES
  138. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  139. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  140. #endif
  141. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  142. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  143. #endif
  144. #endif
  145.  
  146. #ifndef MULTIPLY        /* default definition */
  147. #define MULTIPLY(var,const)  ((var) * (const))
  148. #endif
  149.  
  150. /* Precomputed idct value arrays. */
  151.  
  152. static DCTELEM PreIDCT[64][64];
  153.  
  154. /* Pre compute singleton coefficient IDCT values. */
  155. void
  156. init_pre_idct() {
  157.   int i;
  158.  
  159.   for (i=0; i<64; i++) {
  160.     memset((char *) PreIDCT[i], 0, 64*sizeof(DCTELEM));
  161.     PreIDCT[i][i] = 2048;
  162.     mpeg_j_rev_dct(PreIDCT[i]);
  163.   }
  164. }
  165.  
  166. #ifndef ORIG_DCT
  167.   
  168.  
  169. /*
  170.  * Perform the inverse DCT on one block of coefficients.
  171.  */
  172.  
  173. void
  174. j_rev_dct_sparse (data, pos)
  175.      DCTBLOCK data;
  176.      int pos;
  177. {
  178.   register DCTELEM *dataptr;
  179.   short int val;
  180.   DCTELEM *ndataptr;
  181.   int scale, coeff, rr;
  182.   register int *dp;
  183.   register int v;
  184.  
  185.   /* If DC Coefficient. */
  186.   
  187.   if (pos == 0) {
  188.     dp = (int *)data;
  189.     v = *data;
  190.     /* Compute 32 bit value to assign.  This speeds things up a bit */
  191.     if (v < 0) val = (v-3)>>3;
  192.     else val = (v+4)>>3;
  193.     v = val | (val << 16);
  194.     dp[0] = v;      dp[1] = v;      dp[2] = v;      dp[3] = v;
  195.     dp[4] = v;      dp[5] = v;      dp[6] = v;      dp[7] = v;
  196.     dp[8] = v;      dp[9] = v;      dp[10] = v;      dp[11] = v;
  197.     dp[12] = v;      dp[13] = v;      dp[14] = v;      dp[15] = v;
  198.     dp[16] = v;      dp[17] = v;      dp[18] = v;      dp[19] = v;
  199.     dp[20] = v;      dp[21] = v;      dp[22] = v;      dp[23] = v;
  200.     dp[24] = v;      dp[25] = v;      dp[26] = v;      dp[27] = v;
  201.     dp[28] = v;      dp[29] = v;      dp[30] = v;      dp[31] = v;
  202.     return;
  203.   }
  204.   
  205.   /* Some other coefficient. */
  206.   dataptr = (DCTELEM *)data;
  207.   coeff = dataptr[pos];
  208.   ndataptr = PreIDCT[pos];
  209.  
  210.   for (rr=0; rr<4; rr++) {
  211.     dataptr[0] = (ndataptr[0] * coeff) >> (CONST_BITS-2);
  212.     dataptr[1] = (ndataptr[1] * coeff) >> (CONST_BITS-2);
  213.     dataptr[2] = (ndataptr[2] * coeff) >> (CONST_BITS-2);
  214.     dataptr[3] = (ndataptr[3] * coeff) >> (CONST_BITS-2);
  215.     dataptr[4] = (ndataptr[4] * coeff) >> (CONST_BITS-2);
  216.     dataptr[5] = (ndataptr[5] * coeff) >> (CONST_BITS-2);
  217.     dataptr[6] = (ndataptr[6] * coeff) >> (CONST_BITS-2);
  218.     dataptr[7] = (ndataptr[7] * coeff) >> (CONST_BITS-2);
  219.     dataptr[8] = (ndataptr[8] * coeff) >> (CONST_BITS-2);
  220.     dataptr[9] = (ndataptr[9] * coeff) >> (CONST_BITS-2);
  221.     dataptr[10] = (ndataptr[10] * coeff) >> (CONST_BITS-2);
  222.     dataptr[11] = (ndataptr[11] * coeff) >> (CONST_BITS-2);
  223.     dataptr[12] = (ndataptr[12] * coeff) >> (CONST_BITS-2);
  224.     dataptr[13] = (ndataptr[13] * coeff) >> (CONST_BITS-2);
  225.     dataptr[14] = (ndataptr[14] * coeff) >> (CONST_BITS-2);
  226.     dataptr[15] = (ndataptr[15] * coeff) >> (CONST_BITS-2);
  227.     dataptr += 16;
  228.     ndataptr += 16;
  229.   }
  230.   return;
  231. }
  232.  
  233.  
  234. void
  235. mpeg_j_rev_dct (data)
  236.      DCTBLOCK data;
  237. {
  238.   INT32 tmp0, tmp1, tmp2, tmp3;
  239.   INT32 tmp10, tmp11, tmp12, tmp13;
  240.   INT32 z1, z2, z3, z4, z5;
  241.   INT32 d0, d1, d2, d3, d4, d5, d6, d7;
  242.   register DCTELEM *dataptr;
  243.   int rowctr;
  244.   SHIFT_TEMPS
  245.    
  246.   /* Pass 1: process rows. */
  247.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  248.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  249.  
  250.   dataptr = data;
  251.  
  252.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  253.     /* Due to quantization, we will usually find that many of the input
  254.      * coefficients are zero, especially the AC terms.  We can exploit this
  255.      * by short-circuiting the IDCT calculation for any row in which all
  256.      * the AC terms are zero.  In that case each output is equal to the
  257.      * DC coefficient (with scale factor as needed).
  258.      * With typical images and quantization tables, half or more of the
  259.      * row DCT calculations can be simplified this way.
  260.      */
  261.  
  262.     register int *idataptr = (int*)dataptr;
  263.     d0 = dataptr[0];
  264.     d1 = dataptr[1];
  265.     if ((d1 == 0) && (idataptr[1] | idataptr[2] | idataptr[3]) == 0) {
  266.       /* AC terms all zero */
  267.       if (d0) {
  268.       /* Compute a 32 bit value to assign. */
  269.       DCTELEM dcval = (DCTELEM) (d0 << PASS1_BITS);
  270.       register int v = (dcval & 0xffff) | ((dcval << 16) & 0xffff0000);
  271.       
  272.       idataptr[0] = v;
  273.       idataptr[1] = v;
  274.       idataptr[2] = v;
  275.       idataptr[3] = v;
  276.       }
  277.       
  278.       dataptr += DCTSIZE;    /* advance pointer to next row */
  279.       continue;
  280.     }
  281.     d2 = dataptr[2];
  282.     d3 = dataptr[3];
  283.     d4 = dataptr[4];
  284.     d5 = dataptr[5];
  285.     d6 = dataptr[6];
  286.     d7 = dataptr[7];
  287.  
  288.     /* Even part: reverse the even part of the forward DCT. */
  289.     /* The rotator is sqrt(2)*c(-6). */
  290.     if (d6) {
  291.     if (d4) {
  292.         if (d2) {
  293.         if (d0) {
  294.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  295.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  296.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  297.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  298.  
  299.             tmp0 = (d0 + d4) << CONST_BITS;
  300.             tmp1 = (d0 - d4) << CONST_BITS;
  301.  
  302.             tmp10 = tmp0 + tmp3;
  303.             tmp13 = tmp0 - tmp3;
  304.             tmp11 = tmp1 + tmp2;
  305.             tmp12 = tmp1 - tmp2;
  306.         } else {
  307.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  308.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  309.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  310.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  311.  
  312.             tmp0 = d4 << CONST_BITS;
  313.  
  314.             tmp10 = tmp0 + tmp3;
  315.             tmp13 = tmp0 - tmp3;
  316.             tmp11 = tmp2 - tmp0;
  317.             tmp12 = -(tmp0 + tmp2);
  318.         }
  319.         } else {
  320.         if (d0) {
  321.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  322.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  323.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  324.  
  325.             tmp0 = (d0 + d4) << CONST_BITS;
  326.             tmp1 = (d0 - d4) << CONST_BITS;
  327.  
  328.             tmp10 = tmp0 + tmp3;
  329.             tmp13 = tmp0 - tmp3;
  330.             tmp11 = tmp1 + tmp2;
  331.             tmp12 = tmp1 - tmp2;
  332.         } else {
  333.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  334.             tmp2 = MULTIPLY(d6, -FIX(1.306562965));
  335.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  336.  
  337.             tmp0 = d4 << CONST_BITS;
  338.  
  339.             tmp10 = tmp0 + tmp3;
  340.             tmp13 = tmp0 - tmp3;
  341.             tmp11 = tmp2 - tmp0;
  342.             tmp12 = -(tmp0 + tmp2);
  343.         }
  344.         }
  345.     } else {
  346.         if (d2) {
  347.         if (d0) {
  348.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  349.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  350.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  351.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  352.  
  353.             tmp0 = d0 << CONST_BITS;
  354.  
  355.             tmp10 = tmp0 + tmp3;
  356.             tmp13 = tmp0 - tmp3;
  357.             tmp11 = tmp0 + tmp2;
  358.             tmp12 = tmp0 - tmp2;
  359.         } else {
  360.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  361.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  362.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  363.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  364.  
  365.             tmp10 = tmp3;
  366.             tmp13 = -tmp3;
  367.             tmp11 = tmp2;
  368.             tmp12 = -tmp2;
  369.         }
  370.         } else {
  371.         if (d0) {
  372.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  373.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  374.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  375.  
  376.             tmp0 = d0 << CONST_BITS;
  377.  
  378.             tmp10 = tmp0 + tmp3;
  379.             tmp13 = tmp0 - tmp3;
  380.             tmp11 = tmp0 + tmp2;
  381.             tmp12 = tmp0 - tmp2;
  382.         } else {
  383.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  384.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  385.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  386.  
  387.             tmp10 = tmp3;
  388.             tmp13 = -tmp3;
  389.             tmp11 = tmp2;
  390.             tmp12 = -tmp2;
  391.         }
  392.         }
  393.     }
  394.     } else {
  395.     if (d4) {
  396.         if (d2) {
  397.         if (d0) {
  398.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  399.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  400.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  401.  
  402.             tmp0 = (d0 + d4) << CONST_BITS;
  403.             tmp1 = (d0 - d4) << CONST_BITS;
  404.  
  405.             tmp10 = tmp0 + tmp3;
  406.             tmp13 = tmp0 - tmp3;
  407.             tmp11 = tmp1 + tmp2;
  408.             tmp12 = tmp1 - tmp2;
  409.         } else {
  410.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  411.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  412.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  413.  
  414.             tmp0 = d4 << CONST_BITS;
  415.  
  416.             tmp10 = tmp0 + tmp3;
  417.             tmp13 = tmp0 - tmp3;
  418.             tmp11 = tmp2 - tmp0;
  419.             tmp12 = -(tmp0 + tmp2);
  420.         }
  421.         } else {
  422.         if (d0) {
  423.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  424.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  425.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  426.         } else {
  427.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  428.             tmp10 = tmp13 = d4 << CONST_BITS;
  429.             tmp11 = tmp12 = -tmp10;
  430.         }
  431.         }
  432.     } else {
  433.         if (d2) {
  434.         if (d0) {
  435.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  436.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  437.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  438.  
  439.             tmp0 = d0 << CONST_BITS;
  440.  
  441.             tmp10 = tmp0 + tmp3;
  442.             tmp13 = tmp0 - tmp3;
  443.             tmp11 = tmp0 + tmp2;
  444.             tmp12 = tmp0 - tmp2;
  445.         } else {
  446.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  447.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  448.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  449.  
  450.             tmp10 = tmp3;
  451.             tmp13 = -tmp3;
  452.             tmp11 = tmp2;
  453.             tmp12 = -tmp2;
  454.         }
  455.         } else {
  456.         if (d0) {
  457.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  458.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  459.         } else {
  460.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  461.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  462.         }
  463.         }
  464.     }
  465.     }
  466.  
  467.  
  468.     /* Odd part per figure 8; the matrix is unitary and hence its
  469.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  470.      */
  471.  
  472.     if (d7) {
  473.     if (d5) {
  474.         if (d3) {
  475.         if (d1) {
  476.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  477.             z1 = d7 + d1;
  478.             z2 = d5 + d3;
  479.             z3 = d7 + d3;
  480.             z4 = d5 + d1;
  481.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  482.             
  483.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  484.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  485.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  486.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  487.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  488.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  489.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  490.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  491.             
  492.             z3 += z5;
  493.             z4 += z5;
  494.             
  495.             tmp0 += z1 + z3;
  496.             tmp1 += z2 + z4;
  497.             tmp2 += z2 + z3;
  498.             tmp3 += z1 + z4;
  499.         } else {
  500.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  501.             z1 = d7;
  502.             z2 = d5 + d3;
  503.             z3 = d7 + d3;
  504.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  505.             
  506.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  507.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  508.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  509.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  510.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  511.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  512.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  513.             
  514.             z3 += z5;
  515.             z4 += z5;
  516.             
  517.             tmp0 += z1 + z3;
  518.             tmp1 += z2 + z4;
  519.             tmp2 += z2 + z3;
  520.             tmp3 = z1 + z4;
  521.         }
  522.         } else {
  523.         if (d1) {
  524.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  525.             z1 = d7 + d1;
  526.             z2 = d5;
  527.             z3 = d7;
  528.             z4 = d5 + d1;
  529.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  530.             
  531.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  532.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  533.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  534.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  535.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  536.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  537.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  538.             
  539.             z3 += z5;
  540.             z4 += z5;
  541.             
  542.             tmp0 += z1 + z3;
  543.             tmp1 += z2 + z4;
  544.             tmp2 = z2 + z3;
  545.             tmp3 += z1 + z4;
  546.         } else {
  547.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  548.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  549.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  550.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  551.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  552.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  553.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  554.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  555.             
  556.             z3 += z5;
  557.             z4 += z5;
  558.             
  559.             tmp0 += z3;
  560.             tmp1 += z4;
  561.             tmp2 = z2 + z3;
  562.             tmp3 = z1 + z4;
  563.         }
  564.         }
  565.     } else {
  566.         if (d3) {
  567.         if (d1) {
  568.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  569.             z1 = d7 + d1;
  570.             z3 = d7 + d3;
  571.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  572.             
  573.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  574.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  575.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  576.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  577.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  578.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  579.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  580.             
  581.             z3 += z5;
  582.             z4 += z5;
  583.             
  584.             tmp0 += z1 + z3;
  585.             tmp1 = z2 + z4;
  586.             tmp2 += z2 + z3;
  587.             tmp3 += z1 + z4;
  588.         } else {
  589.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  590.             z3 = d7 + d3;
  591.             
  592.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  593.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  594.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  595.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  596.             z5 = MULTIPLY(z3, FIX(1.175875602));
  597.             z3 = MULTIPLY(z3, - FIX(0.785694958));
  598.             
  599.             tmp0 += z3;
  600.             tmp1 = z2 + z5;
  601.             tmp2 += z3;
  602.             tmp3 = z1 + z5;
  603.         }
  604.         } else {
  605.         if (d1) {
  606.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  607.             z1 = d7 + d1;
  608.             z5 = MULTIPLY(z1, FIX(1.175875602));
  609.  
  610.             z1 = MULTIPLY(z1, FIX(0.275899379));
  611.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  612.             tmp0 = MULTIPLY(d7, - FIX(1.662939224)); 
  613.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  614.             tmp3 = MULTIPLY(d1, FIX(1.111140466));
  615.  
  616.             tmp0 += z1;
  617.             tmp1 = z4 + z5;
  618.             tmp2 = z3 + z5;
  619.             tmp3 += z1;
  620.         } else {
  621.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  622.             tmp0 = MULTIPLY(d7, - FIX(1.387039845));
  623.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  624.             tmp2 = MULTIPLY(d7, - FIX(0.785694958));
  625.             tmp3 = MULTIPLY(d7, FIX(0.275899379));
  626.         }
  627.         }
  628.     }
  629.     } else {
  630.     if (d5) {
  631.         if (d3) {
  632.         if (d1) {
  633.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  634.             z2 = d5 + d3;
  635.             z4 = d5 + d1;
  636.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  637.             
  638.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  639.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  640.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  641.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  642.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  643.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  644.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  645.             
  646.             z3 += z5;
  647.             z4 += z5;
  648.             
  649.             tmp0 = z1 + z3;
  650.             tmp1 += z2 + z4;
  651.             tmp2 += z2 + z3;
  652.             tmp3 += z1 + z4;
  653.         } else {
  654.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  655.             z2 = d5 + d3;
  656.             
  657.             z5 = MULTIPLY(z2, FIX(1.175875602));
  658.             tmp1 = MULTIPLY(d5, FIX(1.662939225));
  659.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  660.             z2 = MULTIPLY(z2, - FIX(1.387039845));
  661.             tmp2 = MULTIPLY(d3, FIX(1.111140466));
  662.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  663.             
  664.             tmp0 = z3 + z5;
  665.             tmp1 += z2;
  666.             tmp2 += z2;
  667.             tmp3 = z4 + z5;
  668.         }
  669.         } else {
  670.         if (d1) {
  671.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  672.             z4 = d5 + d1;
  673.             
  674.             z5 = MULTIPLY(z4, FIX(1.175875602));
  675.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  676.             tmp3 = MULTIPLY(d1, FIX(0.601344887));
  677.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  678.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  679.             z4 = MULTIPLY(z4, FIX(0.785694958));
  680.             
  681.             tmp0 = z1 + z5;
  682.             tmp1 += z4;
  683.             tmp2 = z2 + z5;
  684.             tmp3 += z4;
  685.         } else {
  686.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  687.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  688.             tmp1 = MULTIPLY(d5, FIX(0.275899380));
  689.             tmp2 = MULTIPLY(d5, - FIX(1.387039845));
  690.             tmp3 = MULTIPLY(d5, FIX(0.785694958));
  691.         }
  692.         }
  693.     } else {
  694.         if (d3) {
  695.         if (d1) {
  696.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  697.             z5 = d1 + d3;
  698.             tmp3 = MULTIPLY(d1, FIX(0.211164243));
  699.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  700.             z1 = MULTIPLY(d1, FIX(1.061594337));
  701.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  702.             z4 = MULTIPLY(z5, FIX(0.785694958));
  703.             z5 = MULTIPLY(z5, FIX(1.175875602));
  704.             
  705.             tmp0 = z1 - z4;
  706.             tmp1 = z2 + z4;
  707.             tmp2 += z5;
  708.             tmp3 += z5;
  709.         } else {
  710.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  711.             tmp0 = MULTIPLY(d3, - FIX(0.785694958));
  712.             tmp1 = MULTIPLY(d3, - FIX(1.387039845));
  713.             tmp2 = MULTIPLY(d3, - FIX(0.275899379));
  714.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  715.         }
  716.         } else {
  717.         if (d1) {
  718.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  719.             tmp0 = MULTIPLY(d1, FIX(0.275899379));
  720.             tmp1 = MULTIPLY(d1, FIX(0.785694958));
  721.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  722.             tmp3 = MULTIPLY(d1, FIX(1.387039845));
  723.         } else {
  724.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  725.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  726.         }
  727.         }
  728.     }
  729.     }
  730.  
  731.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  732.  
  733.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  734.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  735.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  736.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  737.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  738.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  739.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  740.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  741.  
  742.     dataptr += DCTSIZE;        /* advance pointer to next row */
  743.   }
  744.  
  745.   /* Pass 2: process columns. */
  746.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  747.   /* and also undo the PASS1_BITS scaling. */
  748.  
  749.   dataptr = data;
  750.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  751.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  752.      * However, the row calculation has created many nonzero AC terms, so the
  753.      * simplification applies less often (typically 5% to 10% of the time).
  754.      * On machines with very fast multiplication, it's possible that the
  755.      * test takes more time than it's worth.  In that case this section
  756.      * may be commented out.
  757.      */
  758.  
  759.     d0 = dataptr[DCTSIZE*0];
  760.     d1 = dataptr[DCTSIZE*1];
  761.     d2 = dataptr[DCTSIZE*2];
  762.     d3 = dataptr[DCTSIZE*3];
  763.     d4 = dataptr[DCTSIZE*4];
  764.     d5 = dataptr[DCTSIZE*5];
  765.     d6 = dataptr[DCTSIZE*6];
  766.     d7 = dataptr[DCTSIZE*7];
  767.  
  768.     /* Even part: reverse the even part of the forward DCT. */
  769.     /* The rotator is sqrt(2)*c(-6). */
  770.     if (d6) {
  771.     if (d4) {
  772.         if (d2) {
  773.         if (d0) {
  774.             /* d0 != 0, d2 != 0, d4 != 0, d6 != 0 */
  775.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  776.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  777.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  778.  
  779.             tmp0 = (d0 + d4) << CONST_BITS;
  780.             tmp1 = (d0 - d4) << CONST_BITS;
  781.  
  782.             tmp10 = tmp0 + tmp3;
  783.             tmp13 = tmp0 - tmp3;
  784.             tmp11 = tmp1 + tmp2;
  785.             tmp12 = tmp1 - tmp2;
  786.         } else {
  787.             /* d0 == 0, d2 != 0, d4 != 0, d6 != 0 */
  788.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  789.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  790.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  791.  
  792.             tmp0 = d4 << CONST_BITS;
  793.  
  794.             tmp10 = tmp0 + tmp3;
  795.             tmp13 = tmp0 - tmp3;
  796.             tmp11 = tmp2 - tmp0;
  797.             tmp12 = -(tmp0 + tmp2);
  798.         }
  799.         } else {
  800.         if (d0) {
  801.             /* d0 != 0, d2 == 0, d4 != 0, d6 != 0 */
  802.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  803.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  804.  
  805.             tmp0 = (d0 + d4) << CONST_BITS;
  806.             tmp1 = (d0 - d4) << CONST_BITS;
  807.  
  808.             tmp10 = tmp0 + tmp3;
  809.             tmp13 = tmp0 - tmp3;
  810.             tmp11 = tmp1 + tmp2;
  811.             tmp12 = tmp1 - tmp2;
  812.         } else {
  813.             /* d0 == 0, d2 == 0, d4 != 0, d6 != 0 */
  814.             tmp2 = MULTIPLY(d6, -FIX(1.306562965));
  815.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  816.  
  817.             tmp0 = d4 << CONST_BITS;
  818.  
  819.             tmp10 = tmp0 + tmp3;
  820.             tmp13 = tmp0 - tmp3;
  821.             tmp11 = tmp2 - tmp0;
  822.             tmp12 = -(tmp0 + tmp2);
  823.         }
  824.         }
  825.     } else {
  826.         if (d2) {
  827.         if (d0) {
  828.             /* d0 != 0, d2 != 0, d4 == 0, d6 != 0 */
  829.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  830.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  831.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  832.  
  833.             tmp0 = d0 << CONST_BITS;
  834.  
  835.             tmp10 = tmp0 + tmp3;
  836.             tmp13 = tmp0 - tmp3;
  837.             tmp11 = tmp0 + tmp2;
  838.             tmp12 = tmp0 - tmp2;
  839.         } else {
  840.             /* d0 == 0, d2 != 0, d4 == 0, d6 != 0 */
  841.             z1 = MULTIPLY(d2 + d6, FIX(0.541196100));
  842.             tmp2 = z1 + MULTIPLY(d6, - FIX(1.847759065));
  843.             tmp3 = z1 + MULTIPLY(d2, FIX(0.765366865));
  844.  
  845.             tmp10 = tmp3;
  846.             tmp13 = -tmp3;
  847.             tmp11 = tmp2;
  848.             tmp12 = -tmp2;
  849.         }
  850.         } else {
  851.         if (d0) {
  852.             /* d0 != 0, d2 == 0, d4 == 0, d6 != 0 */
  853.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  854.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  855.  
  856.             tmp0 = d0 << CONST_BITS;
  857.  
  858.             tmp10 = tmp0 + tmp3;
  859.             tmp13 = tmp0 - tmp3;
  860.             tmp11 = tmp0 + tmp2;
  861.             tmp12 = tmp0 - tmp2;
  862.         } else {
  863.             /* d0 == 0, d2 == 0, d4 == 0, d6 != 0 */
  864.             tmp2 = MULTIPLY(d6, - FIX(1.306562965));
  865.             tmp3 = MULTIPLY(d6, FIX(0.541196100));
  866.  
  867.             tmp10 = tmp3;
  868.             tmp13 = -tmp3;
  869.             tmp11 = tmp2;
  870.             tmp12 = -tmp2;
  871.         }
  872.         }
  873.     }
  874.     } else {
  875.     if (d4) {
  876.         if (d2) {
  877.         if (d0) {
  878.             /* d0 != 0, d2 != 0, d4 != 0, d6 == 0 */
  879.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  880.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  881.  
  882.             tmp0 = (d0 + d4) << CONST_BITS;
  883.             tmp1 = (d0 - d4) << CONST_BITS;
  884.  
  885.             tmp10 = tmp0 + tmp3;
  886.             tmp13 = tmp0 - tmp3;
  887.             tmp11 = tmp1 + tmp2;
  888.             tmp12 = tmp1 - tmp2;
  889.         } else {
  890.             /* d0 == 0, d2 != 0, d4 != 0, d6 == 0 */
  891.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  892.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  893.  
  894.             tmp0 = d4 << CONST_BITS;
  895.  
  896.             tmp10 = tmp0 + tmp3;
  897.             tmp13 = tmp0 - tmp3;
  898.             tmp11 = tmp2 - tmp0;
  899.             tmp12 = -(tmp0 + tmp2);
  900.         }
  901.         } else {
  902.         if (d0) {
  903.             /* d0 != 0, d2 == 0, d4 != 0, d6 == 0 */
  904.             tmp10 = tmp13 = (d0 + d4) << CONST_BITS;
  905.             tmp11 = tmp12 = (d0 - d4) << CONST_BITS;
  906.         } else {
  907.             /* d0 == 0, d2 == 0, d4 != 0, d6 == 0 */
  908.             tmp10 = tmp13 = d4 << CONST_BITS;
  909.             tmp11 = tmp12 = -tmp10;
  910.         }
  911.         }
  912.     } else {
  913.         if (d2) {
  914.         if (d0) {
  915.             /* d0 != 0, d2 != 0, d4 == 0, d6 == 0 */
  916.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  917.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  918.  
  919.             tmp0 = d0 << CONST_BITS;
  920.  
  921.             tmp10 = tmp0 + tmp3;
  922.             tmp13 = tmp0 - tmp3;
  923.             tmp11 = tmp0 + tmp2;
  924.             tmp12 = tmp0 - tmp2;
  925.         } else {
  926.             /* d0 == 0, d2 != 0, d4 == 0, d6 == 0 */
  927.             tmp2 = MULTIPLY(d2, FIX(0.541196100));
  928.             tmp3 = MULTIPLY(d2, FIX(1.306562965));
  929.  
  930.             tmp10 = tmp3;
  931.             tmp13 = -tmp3;
  932.             tmp11 = tmp2;
  933.             tmp12 = -tmp2;
  934.         }
  935.         } else {
  936.         if (d0) {
  937.             /* d0 != 0, d2 == 0, d4 == 0, d6 == 0 */
  938.             tmp10 = tmp13 = tmp11 = tmp12 = d0 << CONST_BITS;
  939.         } else {
  940.             /* d0 == 0, d2 == 0, d4 == 0, d6 == 0 */
  941.             tmp10 = tmp13 = tmp11 = tmp12 = 0;
  942.         }
  943.         }
  944.     }
  945.     }
  946.  
  947.     /* Odd part per figure 8; the matrix is unitary and hence its
  948.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  949.      */
  950.     if (d7) {
  951.     if (d5) {
  952.         if (d3) {
  953.         if (d1) {
  954.             /* d1 != 0, d3 != 0, d5 != 0, d7 != 0 */
  955.             z1 = d7 + d1;
  956.             z2 = d5 + d3;
  957.             z3 = d7 + d3;
  958.             z4 = d5 + d1;
  959.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  960.             
  961.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  962.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  963.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  964.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  965.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  966.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  967.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  968.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  969.             
  970.             z3 += z5;
  971.             z4 += z5;
  972.             
  973.             tmp0 += z1 + z3;
  974.             tmp1 += z2 + z4;
  975.             tmp2 += z2 + z3;
  976.             tmp3 += z1 + z4;
  977.         } else {
  978.             /* d1 == 0, d3 != 0, d5 != 0, d7 != 0 */
  979.             z1 = d7;
  980.             z2 = d5 + d3;
  981.             z3 = d7 + d3;
  982.             z5 = MULTIPLY(z3 + d5, FIX(1.175875602));
  983.             
  984.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  985.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  986.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  987.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  988.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  989.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  990.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  991.             
  992.             z3 += z5;
  993.             z4 += z5;
  994.             
  995.             tmp0 += z1 + z3;
  996.             tmp1 += z2 + z4;
  997.             tmp2 += z2 + z3;
  998.             tmp3 = z1 + z4;
  999.         }
  1000.         } else {
  1001.         if (d1) {
  1002.             /* d1 != 0, d3 == 0, d5 != 0, d7 != 0 */
  1003.             z1 = d7 + d1;
  1004.             z2 = d5;
  1005.             z3 = d7;
  1006.             z4 = d5 + d1;
  1007.             z5 = MULTIPLY(z3 + z4, FIX(1.175875602));
  1008.             
  1009.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1010.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1011.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1012.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1013.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1014.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1015.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1016.             
  1017.             z3 += z5;
  1018.             z4 += z5;
  1019.             
  1020.             tmp0 += z1 + z3;
  1021.             tmp1 += z2 + z4;
  1022.             tmp2 = z2 + z3;
  1023.             tmp3 += z1 + z4;
  1024.         } else {
  1025.             /* d1 == 0, d3 == 0, d5 != 0, d7 != 0 */
  1026.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  1027.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1028.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1029.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  1030.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1031.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1032.             z5 = MULTIPLY(d5 + d7, FIX(1.175875602));
  1033.             
  1034.             z3 += z5;
  1035.             z4 += z5;
  1036.             
  1037.             tmp0 += z3;
  1038.             tmp1 += z4;
  1039.             tmp2 = z2 + z3;
  1040.             tmp3 = z1 + z4;
  1041.         }
  1042.         }
  1043.     } else {
  1044.         if (d3) {
  1045.         if (d1) {
  1046.             /* d1 != 0, d3 != 0, d5 == 0, d7 != 0 */
  1047.             z1 = d7 + d1;
  1048.             z3 = d7 + d3;
  1049.             z5 = MULTIPLY(z3 + d1, FIX(1.175875602));
  1050.             
  1051.             tmp0 = MULTIPLY(d7, FIX(0.298631336)); 
  1052.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1053.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1054.             z1 = MULTIPLY(z1, - FIX(0.899976223));
  1055.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1056.             z3 = MULTIPLY(z3, - FIX(1.961570560));
  1057.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1058.             
  1059.             z3 += z5;
  1060.             z4 += z5;
  1061.             
  1062.             tmp0 += z1 + z3;
  1063.             tmp1 = z2 + z4;
  1064.             tmp2 += z2 + z3;
  1065.             tmp3 += z1 + z4;
  1066.         } else {
  1067.             /* d1 == 0, d3 != 0, d5 == 0, d7 != 0 */
  1068.             z3 = d7 + d3;
  1069.             
  1070.             tmp0 = MULTIPLY(d7, - FIX(0.601344887)); 
  1071.             z1 = MULTIPLY(d7, - FIX(0.899976223));
  1072.             tmp2 = MULTIPLY(d3, FIX(0.509795579));
  1073.             z2 = MULTIPLY(d3, - FIX(2.562915447));
  1074.             z5 = MULTIPLY(z3, FIX(1.175875602));
  1075.             z3 = MULTIPLY(z3, - FIX(0.785694958));
  1076.             
  1077.             tmp0 += z3;
  1078.             tmp1 = z2 + z5;
  1079.             tmp2 += z3;
  1080.             tmp3 = z1 + z5;
  1081.         }
  1082.         } else {
  1083.         if (d1) {
  1084.             /* d1 != 0, d3 == 0, d5 == 0, d7 != 0 */
  1085.             z1 = d7 + d1;
  1086.             z5 = MULTIPLY(z1, FIX(1.175875602));
  1087.  
  1088.             z1 = MULTIPLY(z1, FIX(0.275899379));
  1089.             z3 = MULTIPLY(d7, - FIX(1.961570560));
  1090.             tmp0 = MULTIPLY(d7, - FIX(1.662939224)); 
  1091.             z4 = MULTIPLY(d1, - FIX(0.390180644));
  1092.             tmp3 = MULTIPLY(d1, FIX(1.111140466));
  1093.  
  1094.             tmp0 += z1;
  1095.             tmp1 = z4 + z5;
  1096.             tmp2 = z3 + z5;
  1097.             tmp3 += z1;
  1098.         } else {
  1099.             /* d1 == 0, d3 == 0, d5 == 0, d7 != 0 */
  1100.             tmp0 = MULTIPLY(d7, - FIX(1.387039845));
  1101.             tmp1 = MULTIPLY(d7, FIX(1.175875602));
  1102.             tmp2 = MULTIPLY(d7, - FIX(0.785694958));
  1103.             tmp3 = MULTIPLY(d7, FIX(0.275899379));
  1104.         }
  1105.         }
  1106.     }
  1107.     } else {
  1108.     if (d5) {
  1109.         if (d3) {
  1110.         if (d1) {
  1111.             /* d1 != 0, d3 != 0, d5 != 0, d7 == 0 */
  1112.             z2 = d5 + d3;
  1113.             z4 = d5 + d1;
  1114.             z5 = MULTIPLY(d3 + z4, FIX(1.175875602));
  1115.             
  1116.             tmp1 = MULTIPLY(d5, FIX(2.053119869));
  1117.             tmp2 = MULTIPLY(d3, FIX(3.072711026));
  1118.             tmp3 = MULTIPLY(d1, FIX(1.501321110));
  1119.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1120.             z2 = MULTIPLY(z2, - FIX(2.562915447));
  1121.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1122.             z4 = MULTIPLY(z4, - FIX(0.390180644));
  1123.             
  1124.             z3 += z5;
  1125.             z4 += z5;
  1126.             
  1127.             tmp0 = z1 + z3;
  1128.             tmp1 += z2 + z4;
  1129.             tmp2 += z2 + z3;
  1130.             tmp3 += z1 + z4;
  1131.         } else {
  1132.             /* d1 == 0, d3 != 0, d5 != 0, d7 == 0 */
  1133.             z2 = d5 + d3;
  1134.             
  1135.             z5 = MULTIPLY(z2, FIX(1.175875602));
  1136.             tmp1 = MULTIPLY(d5, FIX(1.662939225));
  1137.             z4 = MULTIPLY(d5, - FIX(0.390180644));
  1138.             z2 = MULTIPLY(z2, - FIX(1.387039845));
  1139.             tmp2 = MULTIPLY(d3, FIX(1.111140466));
  1140.             z3 = MULTIPLY(d3, - FIX(1.961570560));
  1141.             
  1142.             tmp0 = z3 + z5;
  1143.             tmp1 += z2;
  1144.             tmp2 += z2;
  1145.             tmp3 = z4 + z5;
  1146.         }
  1147.         } else {
  1148.         if (d1) {
  1149.             /* d1 != 0, d3 == 0, d5 != 0, d7 == 0 */
  1150.             z4 = d5 + d1;
  1151.             
  1152.             z5 = MULTIPLY(z4, FIX(1.175875602));
  1153.             z1 = MULTIPLY(d1, - FIX(0.899976223));
  1154.             tmp3 = MULTIPLY(d1, FIX(0.601344887));
  1155.             tmp1 = MULTIPLY(d5, - FIX(0.509795578));
  1156.             z2 = MULTIPLY(d5, - FIX(2.562915447));
  1157.             z4 = MULTIPLY(z4, FIX(0.785694958));
  1158.             
  1159.             tmp0 = z1 + z5;
  1160.             tmp1 += z4;
  1161.             tmp2 = z2 + z5;
  1162.             tmp3 += z4;
  1163.         } else {
  1164.             /* d1 == 0, d3 == 0, d5 != 0, d7 == 0 */
  1165.             tmp0 = MULTIPLY(d5, FIX(1.175875602));
  1166.             tmp1 = MULTIPLY(d5, FIX(0.275899380));
  1167.             tmp2 = MULTIPLY(d5, - FIX(1.387039845));
  1168.             tmp3 = MULTIPLY(d5, FIX(0.785694958));
  1169.         }
  1170.         }
  1171.     } else {
  1172.         if (d3) {
  1173.         if (d1) {
  1174.             /* d1 != 0, d3 != 0, d5 == 0, d7 == 0 */
  1175.             z5 = d1 + d3;
  1176.             tmp3 = MULTIPLY(d1, FIX(0.211164243));
  1177.             tmp2 = MULTIPLY(d3, - FIX(1.451774981));
  1178.             z1 = MULTIPLY(d1, FIX(1.061594337));
  1179.             z2 = MULTIPLY(d3, - FIX(2.172734803));
  1180.             z4 = MULTIPLY(z5, FIX(0.785694958));
  1181.             z5 = MULTIPLY(z5, FIX(1.175875602));
  1182.             
  1183.             tmp0 = z1 - z4;
  1184.             tmp1 = z2 + z4;
  1185.             tmp2 += z5;
  1186.             tmp3 += z5;
  1187.         } else {
  1188.             /* d1 == 0, d3 != 0, d5 == 0, d7 == 0 */
  1189.             tmp0 = MULTIPLY(d3, - FIX(0.785694958));
  1190.             tmp1 = MULTIPLY(d3, - FIX(1.387039845));
  1191.             tmp2 = MULTIPLY(d3, - FIX(0.275899379));
  1192.             tmp3 = MULTIPLY(d3, FIX(1.175875602));
  1193.         }
  1194.         } else {
  1195.         if (d1) {
  1196.             /* d1 != 0, d3 == 0, d5 == 0, d7 == 0 */
  1197.             tmp0 = MULTIPLY(d1, FIX(0.275899379));
  1198.             tmp1 = MULTIPLY(d1, FIX(0.785694958));
  1199.             tmp2 = MULTIPLY(d1, FIX(1.175875602));
  1200.             tmp3 = MULTIPLY(d1, FIX(1.387039845));
  1201.         } else {
  1202.             /* d1 == 0, d3 == 0, d5 == 0, d7 == 0 */
  1203.             tmp0 = tmp1 = tmp2 = tmp3 = 0;
  1204.         }
  1205.         }
  1206.     }
  1207.     }
  1208.  
  1209.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1210.  
  1211.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1212.                        CONST_BITS+PASS1_BITS+3);
  1213.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1214.                        CONST_BITS+PASS1_BITS+3);
  1215.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1216.                        CONST_BITS+PASS1_BITS+3);
  1217.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1218.                        CONST_BITS+PASS1_BITS+3);
  1219.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1220.                        CONST_BITS+PASS1_BITS+3);
  1221.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1222.                        CONST_BITS+PASS1_BITS+3);
  1223.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1224.                        CONST_BITS+PASS1_BITS+3);
  1225.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1226.                        CONST_BITS+PASS1_BITS+3);
  1227.     
  1228.     dataptr++;            /* advance pointer to next column */
  1229.   }
  1230. }
  1231.  
  1232. #else
  1233.  
  1234.  
  1235. void
  1236. j_rev_dct_sparse (data, pos) 
  1237.      DCTBLOCK data;
  1238.      int pos;
  1239. {
  1240.   j_rev_dct(data);
  1241. }
  1242.  
  1243. void
  1244. j_rev_dct (data)
  1245.   DCTBLOCK data;
  1246. {
  1247.   INT32 tmp0, tmp1, tmp2, tmp3;
  1248.   INT32 tmp10, tmp11, tmp12, tmp13;
  1249.   INT32 z1, z2, z3, z4, z5;
  1250.   register DCTELEM *dataptr;
  1251.   int rowctr;
  1252.   SHIFT_TEMPS
  1253.  
  1254.   /* Pass 1: process rows. */
  1255.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  1256.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  1257.  
  1258.   dataptr = data;
  1259.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1260.     /* Due to quantization, we will usually find that many of the input
  1261.      * coefficients are zero, especially the AC terms.  We can exploit this
  1262.      * by short-circuiting the IDCT calculation for any row in which all
  1263.      * the AC terms are zero.  In that case each output is equal to the
  1264.      * DC coefficient (with scale factor as needed).
  1265.      * With typical images and quantization tables, half or more of the
  1266.      * row DCT calculations can be simplified this way.
  1267.      */
  1268.  
  1269.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  1270.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  1271.       /* AC terms all zero */
  1272.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  1273.       
  1274.       dataptr[0] = dcval;
  1275.       dataptr[1] = dcval;
  1276.       dataptr[2] = dcval;
  1277.       dataptr[3] = dcval;
  1278.       dataptr[4] = dcval;
  1279.       dataptr[5] = dcval;
  1280.       dataptr[6] = dcval;
  1281.       dataptr[7] = dcval;
  1282.       
  1283.       dataptr += DCTSIZE;    /* advance pointer to next row */
  1284.       continue;
  1285.     }
  1286.  
  1287.     /* Even part: reverse the even part of the forward DCT. */
  1288.     /* The rotator is sqrt(2)*c(-6). */
  1289.  
  1290.     z2 = (INT32) dataptr[2];
  1291.     z3 = (INT32) dataptr[6];
  1292.  
  1293.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1294.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1295.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1296.  
  1297.     tmp0 = ((INT32) dataptr[0] + (INT32) dataptr[4]) << CONST_BITS;
  1298.     tmp1 = ((INT32) dataptr[0] - (INT32) dataptr[4]) << CONST_BITS;
  1299.  
  1300.     tmp10 = tmp0 + tmp3;
  1301.     tmp13 = tmp0 - tmp3;
  1302.     tmp11 = tmp1 + tmp2;
  1303.     tmp12 = tmp1 - tmp2;
  1304.     
  1305.     /* Odd part per figure 8; the matrix is unitary and hence its
  1306.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1307.      */
  1308.  
  1309.     tmp0 = (INT32) dataptr[7];
  1310.     tmp1 = (INT32) dataptr[5];
  1311.     tmp2 = (INT32) dataptr[3];
  1312.     tmp3 = (INT32) dataptr[1];
  1313.  
  1314.     z1 = tmp0 + tmp3;
  1315.     z2 = tmp1 + tmp2;
  1316.     z3 = tmp0 + tmp2;
  1317.     z4 = tmp1 + tmp3;
  1318.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1319.     
  1320.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1321.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1322.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1323.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1324.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1325.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1326.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1327.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1328.     
  1329.     z3 += z5;
  1330.     z4 += z5;
  1331.     
  1332.     tmp0 += z1 + z3;
  1333.     tmp1 += z2 + z4;
  1334.     tmp2 += z2 + z3;
  1335.     tmp3 += z1 + z4;
  1336.  
  1337.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1338.  
  1339.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  1340.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  1341.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  1342.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  1343.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  1344.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  1345.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  1346.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  1347.  
  1348.     dataptr += DCTSIZE;        /* advance pointer to next row */
  1349.   }
  1350.  
  1351.   /* Pass 2: process columns. */
  1352.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  1353.   /* and also undo the PASS1_BITS scaling. */
  1354.  
  1355.   dataptr = data;
  1356.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  1357.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  1358.      * However, the row calculation has created many nonzero AC terms, so the
  1359.      * simplification applies less often (typically 5% to 10% of the time).
  1360.      * On machines with very fast multiplication, it's possible that the
  1361.      * test takes more time than it's worth.  In that case this section
  1362.      * may be commented out.
  1363.      */
  1364.  
  1365. #ifndef NO_ZERO_COLUMN_TEST
  1366.     if ((dataptr[DCTSIZE*1] | dataptr[DCTSIZE*2] | dataptr[DCTSIZE*3] |
  1367.      dataptr[DCTSIZE*4] | dataptr[DCTSIZE*5] | dataptr[DCTSIZE*6] |
  1368.      dataptr[DCTSIZE*7]) == 0) {
  1369.       /* AC terms all zero */
  1370.       DCTELEM dcval = (DCTELEM) DESCALE((INT32) dataptr[0], PASS1_BITS+3);
  1371.       
  1372.       dataptr[DCTSIZE*0] = dcval;
  1373.       dataptr[DCTSIZE*1] = dcval;
  1374.       dataptr[DCTSIZE*2] = dcval;
  1375.       dataptr[DCTSIZE*3] = dcval;
  1376.       dataptr[DCTSIZE*4] = dcval;
  1377.       dataptr[DCTSIZE*5] = dcval;
  1378.       dataptr[DCTSIZE*6] = dcval;
  1379.       dataptr[DCTSIZE*7] = dcval;
  1380.       
  1381.       dataptr++;        /* advance pointer to next column */
  1382.       continue;
  1383.     }
  1384. #endif
  1385.  
  1386.     /* Even part: reverse the even part of the forward DCT. */
  1387.     /* The rotator is sqrt(2)*c(-6). */
  1388.  
  1389.     z2 = (INT32) dataptr[DCTSIZE*2];
  1390.     z3 = (INT32) dataptr[DCTSIZE*6];
  1391.  
  1392.     z1 = MULTIPLY(z2 + z3, FIX(0.541196100));
  1393.     tmp2 = z1 + MULTIPLY(z3, - FIX(1.847759065));
  1394.     tmp3 = z1 + MULTIPLY(z2, FIX(0.765366865));
  1395.  
  1396.     tmp0 = ((INT32) dataptr[DCTSIZE*0] + (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1397.     tmp1 = ((INT32) dataptr[DCTSIZE*0] - (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  1398.  
  1399.     tmp10 = tmp0 + tmp3;
  1400.     tmp13 = tmp0 - tmp3;
  1401.     tmp11 = tmp1 + tmp2;
  1402.     tmp12 = tmp1 - tmp2;
  1403.     
  1404.     /* Odd part per figure 8; the matrix is unitary and hence its
  1405.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  1406.      */
  1407.  
  1408.     tmp0 = (INT32) dataptr[DCTSIZE*7];
  1409.     tmp1 = (INT32) dataptr[DCTSIZE*5];
  1410.     tmp2 = (INT32) dataptr[DCTSIZE*3];
  1411.     tmp3 = (INT32) dataptr[DCTSIZE*1];
  1412.  
  1413.     z1 = tmp0 + tmp3;
  1414.     z2 = tmp1 + tmp2;
  1415.     z3 = tmp0 + tmp2;
  1416.     z4 = tmp1 + tmp3;
  1417.     z5 = MULTIPLY(z3 + z4, FIX(1.175875602)); /* sqrt(2) * c3 */
  1418.     
  1419.     tmp0 = MULTIPLY(tmp0, FIX(0.298631336)); /* sqrt(2) * (-c1+c3+c5-c7) */
  1420.     tmp1 = MULTIPLY(tmp1, FIX(2.053119869)); /* sqrt(2) * ( c1+c3-c5+c7) */
  1421.     tmp2 = MULTIPLY(tmp2, FIX(3.072711026)); /* sqrt(2) * ( c1+c3+c5-c7) */
  1422.     tmp3 = MULTIPLY(tmp3, FIX(1.501321110)); /* sqrt(2) * ( c1+c3-c5-c7) */
  1423.     z1 = MULTIPLY(z1, - FIX(0.899976223)); /* sqrt(2) * (c7-c3) */
  1424.     z2 = MULTIPLY(z2, - FIX(2.562915447)); /* sqrt(2) * (-c1-c3) */
  1425.     z3 = MULTIPLY(z3, - FIX(1.961570560)); /* sqrt(2) * (-c3-c5) */
  1426.     z4 = MULTIPLY(z4, - FIX(0.390180644)); /* sqrt(2) * (c5-c3) */
  1427.     
  1428.     z3 += z5;
  1429.     z4 += z5;
  1430.     
  1431.     tmp0 += z1 + z3;
  1432.     tmp1 += z2 + z4;
  1433.     tmp2 += z2 + z3;
  1434.     tmp3 += z1 + z4;
  1435.  
  1436.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  1437.  
  1438.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  1439.                        CONST_BITS+PASS1_BITS+3);
  1440.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  1441.                        CONST_BITS+PASS1_BITS+3);
  1442.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  1443.                        CONST_BITS+PASS1_BITS+3);
  1444.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  1445.                        CONST_BITS+PASS1_BITS+3);
  1446.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  1447.                        CONST_BITS+PASS1_BITS+3);
  1448.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  1449.                        CONST_BITS+PASS1_BITS+3);
  1450.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  1451.                        CONST_BITS+PASS1_BITS+3);
  1452.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  1453.                        CONST_BITS+PASS1_BITS+3);
  1454.     
  1455.     dataptr++;            /* advance pointer to next column */
  1456.   }
  1457. }
  1458.  
  1459.  
  1460. #endif
  1461.